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In this work, we have theoreticaUy analyzed and numerically evaluated the accuracy of 
high-order lattice Boltzmann (LB) models for capturing non-equilibrium effects in rar- 
efied gas flows. In the incompressible limit, the LB equation is proved to be equivalent 
to the linearized Bhatnagar-Gross-Krook (BGK) equation. Therefore, when the same 
Gauss-Hermite quadrature is used, LB method closely assembles the discrete velocity 
method (DVM) . In addition, the order of Hermite expansion for the equilibrium distri- 
bution function is found not to be correlated with the approximation order in terms of the 
Knuds en number to the BGK equation, which was previously suggested by 



Shan et al 



(|2006l ). Furthermore, we have numerically evaluated the LB models for a standing-shear- 
wave problem, which is designed specifically for assessing model accuracy by excluding 
the influence of gas molecule/surface interactions at wall boundaries. The numerical sim- 
ulation results confirm that the high-order terms in the discrete equilibrium distribution 
function play a negligible role. Meanwhile, appropriate Gauss-Hermite quadrature has 
the most significant effect on whether LB models can describe the essential flow physics 
of rarefled gas accurately. For the same order of the Gauss-Hermite quadrature, the exact 
abscissae will also modestly influence numerical accuracy. Using the same Gauss-Hermite 
quadrature, the numerical results of both LB and DVM methods are in excellent agree- 
ment for flows across a broad range of the Knudsen numbers, which conflrms that the 
LB simulation is similar to the DVM process. Therefore, LB method can offer flexible 
models suitable for simulating continuum flows at Navier Stokes level and rarefled gas 
flows at the linearized Boltzmann equation level. 
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1. Introduction 

Rarefied gas flows have recently attracted significant research interest due to the rapid 
development of micro/nano-fluidic technologies. Gaseous transport in micro/nano devices 
is often found to be non-equilibri um, and non-equilibrium phenomena have not yet been 



well understood (Ho fc Tai 



19981 ). The conventional theory to describe gas flows is the 



Navier Stokes equations, which assume that the fluid is in a quasi-equilibrium state. 
However, for non-equilibrium flows, the Navier Stokes equations break down because 
that the molecular nature of the gas strongly affects the bulk flow behavior i.e. the 
gas can no longer be regarded as a fluid continuum. Whether gas flows are in local 
equilibrium or not can be classified by the non-dimensional Knudsen number, Kn, defined 
as the ratio of mean free path and the device characteristic length scale. The Navier 
Stokes equations with no- velocity-slip wall boundary condition arc only appropriate when 
Kn < 0.001. However, gas flows in micro/nano-fiuidic devices are often in the slip fiow 
regime (0.001 < Kn < 0.1) or the transition fiow regime (0.1 < Kn < 10). In these 
regimes, the gas flow cannot properly be described as a continuous flow, nor as a free 
molecular flow. In practice, most devices operate with a range of Knudsen numbers in 
different parts of the device; this makes it even more difficult to develop a generalized 
flow model. 

Direct simulation Monte Carlo (DSMC) methods and direct numerical simulation of 
the Boltzmann equation can provide accurate solutions for rarefied gas flows. However, 
these are computationally intractable for 3D fiow systems, and impractical with the 
current computer technology, especially for the low speed gas fiows usually encountered 
in micro/nano-systems. Statistically, one needs to take significantly more samples of the 
flow field at any point for the DSMC method to resolve flows with low Mach numbers. The 
direct simulations based on the Boltzmann equation requires significant computational 
resources for integrating the velocity space ranging from — oo to +oo. In addition, it is 
usually difficult to solve the full Boltzmann equation directly via either numerical or 
analytical methods. 

Meanwhile, the continuum methods beyond the Navier Stokes level have failed to pro- 



duce satisfactory results for gas flows in the transition flow regime (jLockerbv fc Reese 



20081 ). It is well-known that continuum expressions for the viscous stress and heat flux 



in gases may be derived from the fundamental Boltzmann equation via cither a Kn- 
series solution (known as the Chapman-Enskog approach) or by an expansion of the 
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distri bution function as a series of Hermite tensor polynomials (jChapman fc Cowling 



199lh . To first order (i.e. for near-equilibrium flows) both approaches yield the Navier 



Stokes equations. However, the solution methods can be continued to second and higher 
orders, incorporating more and more of the salient characteristics of a non-equilibrium 
flow. The classical second-order stress and heat flux expressions are the Burnett equa- 
tions (from the Chapman-Enskog a pproach), and the Grad 13- moment equations (from 
the Hermite polynomial method) ([Chapman fc Cowling Il99ll ). These can be seen as 
corrections to the Navier Stokes constitutive relations to make them more appropri- 
ate to continuum-transition flows. However, different physical interpretations of the so- 
lution methods at second and higher orders have recently led to a variety of sets o f 



equ ations, including the Bhatn agar-Gross-Krook (B GK)-Burnett(Balakrishnan 



Eu (Al-Ghoul fc ChanEu 



ized moment (R13) 



20041). augmented Bur nett (jZhong et al 



20041) ■ 



19931 ). and regular- 



Struchtrup fc Torrilhonll2003r ) equations. While each purports to be 



the proper high-order correction to the stress and heat flux (there is no disagreement 
about the form of the Navier Stokes e quations at first-order), none of these models 
are satisfactory (jLockerbv fc Reesd 120081 ) . In addition, these models suffer unknown ad- 
ditional boundary conditions at solid walls to appropriately reflect gas molecule/wall 
surface interactions. 

The lattice Boltzmann (LB) approach offers an alternative method for rarefled gas flow 
simulations. Historically, the LB model was evolved fro m the lattice-gas automata (LGA) 
for mimicking the Navi er Stokes hydrodynamics (see 



1998; 



Benzi et al. 



Qian et al. 



1992 : 



Chen fc Doolen 



I992L and references therein). Over the past two decades, the LB 



method has been developed to provide accurate and efficient solutions for continuum 
flow simulations as the validity of the model can be ensured by the Chapman-Enskog 
expansion. Due to its kinetic nature, the LB model has distinct advantages over the con- 
tinuum computational methods, including easy implementation of multi-physical mech- 



models for simulating rarefied gas flows have been demonstrated (e.g. 


Zhang et al. 


2005 


Toschi fc Sued 


2OO5I 


Sbragaglia fc SuccJ 


2OO5I 


20061 


Tang et al. 


20081 


Zhang et al. 


2006 


Shan et al. 


200k 


Ansumali et al. 


2007; 


Kim et al. 


2008 


. 


Yudistiawan et al. 


20081). 



Recently, the LB models were shown to be able to be derived sy stematically from 
the Boltzmann-BGK equation ba sed on the Hermite expansion (see 



Shan fc He 



1998 



He fc Luo 



19976 



Shan et al. 



20061 1. This creates another theoretical foundation different 
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from the Chapman-Enskog expansion, so that higher-order LB approximations to the 
Boltzmann-BGK equation beyond the Navier Stokes level can be constructed by using 
the high-order Hermite expansion with appropriate quadratures. This indicates that high- 
order LB models have the potential to capture non-equilibrium e ffects in rarefied fl ows. 
In addition to the systemic framework of constructing LB models, 



Shan et al. 



( 20061 ) also 



established the link between the orders of Hermite polynomials and Chapman-Enskog 
expansion. The authors concluded that the order of Hermite expansion is responsible for 
obtaining correct velocity moments. The precise relation among the orders of Hermite 
ex pansion, Ch a pman -Enskog expansion and velocity moments was described by Eq.(4.7) 



m 



Shan et al. 



(|2006r ). For instance, the third-order expansion is required for accurate 
pressure tensor and momentum at the Navier Stokes level. These conclusions are key to 
constructing appropriate LB models for non-equilibrium gas flows. However, the numeri- 
cal simulations do not support these conclusions. In contrast, the simulation data showed 
that the higher order terms in the e g uilibrium distrib ution function have negligible in- 



fluence for low speed rarefied flows (jKim et al. 



2008|). This indicates that the Hermite 



expansion order is not related to the o rder of Chapman- Enksog expansion, in contrary 



to the theoretical conclusions drawn by 



Shan et al. 



pood) . 



In this work, we aim to answer this question whether the Hermite expansion order is 
important for the LB method, as it is for the Grad's moment method, to capture non- 
equilibrium effects in rarefied flows, especially at micro/nano-scales. Furthermore, we will 
analyze theoretically and numerically the mechanisms that are important in constructing 
high-order LB mo dels for rarefled ga s dynamics. We will discuss the differences between 
the approaches of lShan et al\ (|2006[ ) and Grad's moment method. To help us to under- 
stand the modeling capability of the LB method for rarefled gas dynamics, we will also 
analyze the similarities and differences between the LB method and the discrete velocity 
method (DVM) of solving the BGK equation. In particular, we will prove that the Her- 
mite expansion order is not important for the flows that the linearized BGK equation can 
accurately describe. Since the important nonlinear constitu tive relations in the Knudsen 
layer are still not captured satisfactorily (jTang et aLll2008f ). our num erical analysis will 
be ba sed on a standing-shear-wave problem specifically designed by iLockerbv fc Reese 
(j2008l ) to exclude the effect of gas molecule/wall interactions, so we can concentrate on 
the model capabilities. 
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2. Lattice Boltzmann simulation of rarefied gas flows 

2.1. Lattice Boltzmann equation 

Although the LB models were ori ginally developed f ro m LGA, the link t o the kinetic 
theory has late been estabhshed by He fc Luo (1997^); Shan &: He (1998); Shan et al 



(|2006[ ). Consequently, the LB approach may be cons idered as a special finite differ- 



ence scheme of solving the Boltzmann-BGK equation (|Luo 



2000) . This theoretical hnk 



indicates that the LB methodology may provide a reasonable approximation to the 
Boltzmann-BGK equation. The central question is how accurate the LB models can 
capture non-equilibrium effects in rarefied gas dynamics. To answer this question, we 
will revis it the derivation p rocess of LB models from the Boltzmann-BGK equation pro- 
posed bv lShan et all (j2006f ) and we will emphasize on the model capability in describing 
rarefied gas flows. 

The original Boltzmann-BGK equation can be written as: 



df 



P 



(2.1) 



dt ' ' ■ fj, 

where / denotes the distribution function, ^ the phase velocity, p the pressure, g the 
body force and /i the gas viscosity. Using the well-known Chapman-Enskog expansion, 
the collision frequency can be represented by the ratio of pressure and gas viscosity, 
which is convenient to obtain the Knudsen number definition consistent with that of 
hydrodynamic models. Without losing generality, we define the following non-dimensional 
variables: 



f ~ Or, u 



,t. 



rh 



i ^ T 



--,T = 



(2.2) 



where u is the macroscopic velocity, R the gas constant, T the gas temperature, To the 
reference temperature, r the spatial position and 6 the inverse of the characteristic length 
of the flow system. The symbol hat, which denotes a dimensionless value, will hereinafter 
be omitted. We define the Knudsen number using macroscopic properties as below: 



Kn ^ 



P 



(2.3) 



By using these non-dimensional variables, the non-dimensional form of the Boltzmann- 
BGK equation becomes 



f +l-v/ + ,.v,/^-^(/-/-). 



(2.4) 
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where the Maxweh distribution in Z3-dimensional Cartesian coordinates can be written 
as 



eq 



■ exp 



l21 



2T 



(2.5) 



(27rr)^/2 

From the non-dimensional format of Eq. (|2.4p . we can clearly see the relationship between 
the relaxation time and the mean free path (i. e. Knudsen numb er), which plays a key 
role in LB simulation of rarefied gas flows fe.g. lZhang et aLll2005() . 



To discretize the velocity space, we project the distribution function onto a functional 
space spanned by the orthogonal Hermite basis: 

N 

(2.6) 



^ 1 



ri=0 



where the nth order Hermite polynomial. The weight function uj{£) is given by 

1 



^(0 



(27r)^/2 



and the coefficients a^") are 



J J Aio) 



(n) 

The coefficient a^q for the equilibrium distribution is 



,(«) - 



(2.7) 



(2.8) 



(2.9) 



where Wa and ^q, a = I,-- - ,d, are the weights and abscissae of a Gauss- Hermite 
quadrature of degree ^ 27V respectively. Herein, the distribution function is approxi- 
mated by the first N Hermite polynomial. Using the derivation relation, the body force 
term F{r,^,t) = g - ^^f can be approximated as 



AT 



1 



^, ^ (2.10) 

As an example, the second order approximation of the equilibrium distribution and 
the body force are: 



(Op|i + C-« + ^[(4 



D) 



(2.11) 



Fir, 1 1) - L^ii)p {g-i+ig-mu-i)-g-u}, (2.12) 

where T should be set to unity for isothermal problems and p is constant for incompress- 
ible problems. 
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An appropriate Gauss-Hermite quadrature, see the Appendix in I Shan et al\ (|2006l ) for 
a hst of quadratures, can be chosen to evaluate the integral to obtain a^"^ Consequently, 
Ea. (|2.4p can be discretized as 



9o 



(2.13) 



where /a = Waf{r,ia,t)/i^{ia), fa'' = Waf'"^{r,^a,t)/u;{ia) and^Q = WaF{r,^a,t)/Lji^a) 
We have obtained the lattice Boltzmann equation, i.e. Eq. (|2.13p . by discretizing Eq. 
in the velocity space. 



2.2. Numerical schemes, Knudsen number and relaxation time 

An appropriate numerical scheme is now required to solve Eq. (j2.13p . If a finite difference 
scheme is chosen, one can obtain the so-called finite difference lattice Boltzmann model. 
In particular, when the first-order upwind finite-difference scheme is chosen, one can 
obtain the standard form of LB model: 

5t 



Ur + ^aSt,t + St)- fir,t) 



Kn 



(/a - 17) + Stga 



(2.14) 



where the relationship between the relaxation time r and the Knudsen number is estab- 
lished naturally i.e. r — kri/St. For continuum flows where the Navier Stokes equations 
are valid, the above first-order scheme can become effectively second-order accurate in 
both space and time by simply replacing the non-dime nsional relaxation time r with 
T — 0.5 (see 



Reider fc Sterlinelll995 



Sterhng fc CheJl996 !). In doing so, the second order 



discretization error can be absorbed into an artificial viscosity. Therefore, this simple but 
accurate scheme has been widely used to simulate flows at the Navier Stokes level. Since 
any LB model intended to simulate rarefied gas dynamics beyond the Navier Stokes level 
needs to recover the Navier Stokes equation at small Knudsen number, i.e. Kn 0, this 
first-order sch eme with correction has been com monly used in LB simulation of rarefied 
gas flows (e.g. 



Nie et al. 



2n02t 



Tang et aZ.1 120081 ) for rarefied gas problems. However, the 



artificial viscosity has only corrected the momentum transfer to the second order, which 
is only appropriate for the Navier-Stokes hydrodynamics. This correction will lead to 
inconsistency for the transfer of the other higher-order moments, which are essential for 
capturing non-equilibrium effects in rarefied gas fiows. Therefore, the dilemma is that 
we need correction on the relaxation time to recover the Navier Stokes hydrodynamics 
appropriately when the Knudsen number is close to zero where the high-order moments 
are not important. Meanwhile, we should not have this correction for the higher-order 
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moments which are more important to rarefied flows. In deed, the simulation will di- 
verge when the Knudsen number is approaching to 0.5 if no relaxation time correction 
is intro duced. The reaso n is that it goes beyond the stability regime of the relaxation 
( 20021 ) suggested to use the correction when the Knudsen number is 



scheme. 



Lim et 



less than 0.5 and switch to no correction when the Knudsen number is larger than 0.5. 
However, it will lead to inconsistency at the Knudsen number around 0.5 which is the 
most important flow regime in micro/nano-fluidic devices. The above first-order upwind 
scheme should not be used for simulating the gas flows with finite Knudsen numbers. 

To resolve this problem inherited from the standard LB method, we should not use 
the artificial viscosity to achieve correct physics at the Navier Stokes level. We propose 
to discretize Eg. (12. 13 1 using a numerical scheme with second-order accuracy, which was 
(1993) for thermal fiow simulation at the Navier Stokes level: 
5t 



first used by 



He et 



Uir + ^cSt, t + St)- fir, t) = -^ [Uir + ^^St, t + St) ~ /^«(r + ^^6t, t + St)] 



St 

"2kn 
St 



[Ur,t)~ra''{r,t)] 



+ j [9a {r + ^c.St, t + St)+ g^{r, t)] 



By introducing 



fa — fa 

the above implicit scheme can be written as 

St 



—(f - r*)--o 



fair + ^c.St,t + St)-fc.ir,t) = - 



with 



Kn + O.bSt 



/o(r,t)-/^nr,i) 



KnggSt 
Kn + 0.5St 



(2.15) 



(2.16) 



(2.17) 



(2.18) 



PU = ^^afa 



pgSt 



(2.19) 



Therefore, the viscosity is now tRT rather than {t~Q.5)RT. Most importantly, the same 
relation between the relaxation time and the mean free path can be used for the transfer 
of any order moments. 



2.3. High- order lattice Boltzmann models 

Although the construction of LB models based on the Hermite polynomials is straight- 
forward, the Hermite polynomials higher than the third order give irrational roots. The 
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integer stream velocity is an essential feature of LB models, i.e. the simple and efficient 
"stream-collision" mechanism. So high-order LB models, which have non -integer discret e 



velocities, will need additional effort, such as point-wise interpolation (jHe et 



mm- 



Therefore, they essentially become ofF-lattice discrete velocity method for solving the ki- 
netic Boltzmann equation, whic h will increase the computational cost dramatically and 

(j2006l ) suggested a method for searching ab- 



introduce extra numerical error. 



Shan et al. 



scissae on the grid points of Cartesian coordinates to construct high-order LB models 
with intege r discr ete v elocities. The exa mples are D2Qf7 and D2Q21 models given by 



Shan et al. 



mm and 



Kim et 



( 20081 ) (note, we f ollow the conventional terminology 



for the LB models as first introduced bv lQian et al\ (|l992f ) dub bed as DnQm model i.e. 



Chikatamarla fc Karlin 



n dim ensional model with ni discrete velocities). Furthermore 
( 20061 ) proposed an alternative method to seek rational- number approximation to the 
rations of the Hermite roots based on the relation between the entropy and the roots 
of Hermite polynomials. They also proposed the higher-order LB models with integer 
discrete velocity, such as D2Q16 and D2Q25 models. The above high-order LB models 
with integer stream velocities will be numerically examined in this work and the details 
are listed in Table. ([T|). 

Based on the above model construction procedure, the accuracy of LB models depends 
on three level of approximations. Firstly, it depends on the accuracy of the numerical 
scheme for solving Eq. ()2.13p . As we have demonstrated, the commonly used first-order 
upwind scheme will lead to incorrect physics for rarefied flows. Our second-order numer- 
ical scheme given by Eq. (|2.f 7p is essential to capture non-equilibrium effects accurately. 
Secondly, the order o f the Hermite exp ansion was considered to be important to obtain 
the correct moments ( Shan et a/.ll2006 ). Thirdly, the Gauss-Hermite quadrature accuracy 
should be sufficiently high so that the integration of Eq. p.Sp can be evaluated accurately. 
Therefore, the term higher-order LB models here refer to the LB models with high-order 
of Hermite expansion and Gauss-Hermite quadrature in comparison with the standard 
LB model. 



3. Lattice Boltzmann, moment and discrete velocity methods 

3.f . Comparison of Grad's moment method and lattice Boltzmann method 

Similar to the Grad's method for deriving higher order continuum systems (e.g, Grad 13- 
moment equations), using the Hermite expansion to approximate the Boltzmann-BGK 
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Quadrature k 



0,0) 4/9 
V3,0)fs 1/9 
±V3,±^/3) 1/36 



±m, ±m) Wlrn 

±n, ±n) Wl„ 

±m, ±n) W±mW±n 

±n, ±m) W±mW±n 



m = l,n 

T J7 " 

W±n = - 



12(m-^ 



-lOn-'m-'+n* 



12(mii-nii) 



To = (m^ + + - lOn^m^ +n4) /6 



1 


(0,0) 


(575 + 193yT93)/8100 


4 


('■,0)fs 


(3355 - 91^/T93)/18000 


4 


(±r, ±r) 


(655 + 17VT93)/27000 


4 


(±2r, ±2r) 


(685 - 49^/T93)/54000 


4 


(3r, 0)fs 


(1445 - 10lVT93)/162000 



= (125 + 5\/T93)/72 



1 


(0,0) 


91/324 


4 


('■,0)fs 


1/12 


4 


(±r, ±r) 


2/27 


4 


(±2r,0)Fs 


7/360 


4 


(±2r, ±2r) 


1/432 


4 


(3r, 0)fs 


1/1620 



= 3/2 



1 


(0,0) 




4 


(m, 0)fs 




4 


(n, 0)fs 




4 


(±m, ±m) 




4 


(±n, ±n) 




4 


(±m, ±n) 




4 


(±ri, ±m) 





3, n = 7 



9m — 6r7 



- 75m^n^ 

300m^(m^-n^) 
-6m^-27n^m^ + (3n^-2m^)D5 
300n^(n^~m^) 



To = (3m^ + 3n^ +D5)/30 
Ds = \/9m4 - 42n2m2 + 9n'^ 



Table 1. The quadratures of five LB models where k is the number of discrete velocities with 
the same velocity magnitude, the subscript FS denotes a fully symmetric set of points, and 
Woi are the weights. The quadrature accuracy is fifth-order for the D2Q9 model, seventh-order 



of D2Q17 and D2Q21 models can be found in 



1 in 


Shan et al. 


( 


2006 


); 


Kim et al. 


Chikatamarla & Karlin 


( 


2006 


)• 



( 20081 1 while the 
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equation can lead to the LB equation, i.e. Eq. (|2.13|) . However, the major difference is 
that LB models are always staying at the kinetic level, i.e. solving the kinetic equation 
- Eq. (|2.13p . while the Grad's method will produce a set of continuum equations. The 
basic idea of Grad's method is to use the truncated Hermite polynomials to approximate 
the full Boltzmann (or Boltzmann-BGK) equation. Due to the unique feature of Hermite 
polynomial, the moments of up to the chosen truncation order can be described accurately 
by the derived macroscopic moments systems. In contrary, the only explicit effect of the 
truncation on the LB models is on the approximation of the equilibrium distribution 
function and the body force, while the Grad's moment equations do not approximate the 
equilibrium distribution function. 

Although the order of Hermite expansion determines the accuracy level of the moment 
model, which is not the same for the LB models. Essentially, the LB equation i.e. Ea. (|2.13p 
is similar to any model equation which is to simplify the full Boltzmann equation. The 
kinetic process, i.e. gas molecules relaxing to the equilibrium state through collisions, is 
still the same. Therefore, the LB method is very close to the discrete velocity method 
solving the Boltzmann-BGK equation (especially the linearized-BGK equation), which 
we will discuss in the section below. 



3.2. Discrete velocity methods and lattice Boltzmann method 

The above procedure of establishing LB models is similar to the problem solving pro- 
cess of the discrete velocity method, which directly solves the Boltzmann-BGK equa- 
tion. Since DVM has been 



gas dynamics (see 


Mieussens 


2001 


2000 a 


b: Yane & Huane 


1995; 


Aoki et al. 


2002 


1991; 


Valougeorgig 


1988 


; iNaris & Valougeorgis 


2OO5I 


Naris et al. 


200J 




Sharioov 


& Bertoldo 


2009; 


ShariDOv & Kalemoa 


2008, 


and references therein), it is helpful to compare two 



numerical methods in depth. 

The discrete velocity method is to discretize the velocity space bas ed on quadratures 
e.g. Gauss-Hermite quadrature and Newton-Cotes quadrature (see 



Naris et al 



2005 



Naris fc Valougeorgis 



2005 



Valougeorgis 



1988 : 



Yang fc Huang 



19951 ). The first step is 



to non-dimcnsionalize the Boltzmann-BGK equation and obtain the reduced functions, 
e.g. Ga and Qf, in Eqs. (|4.5p and (|4.6p . which are important to reduce computational costs. 
The second step is to apply an appropriate discretization method for the velocity space, 
which is important but difficult because the velocity space ranges from —00 to +00 and 
the properties of conservation and dissipation of the entropy should be kept. A typical 
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choice is the Gauss-Hermite quadrature, which is to be adopted in our simulations. In 
order to reduce the velocity components which need to be integrated from ~oo to +cx), 
curvilinear coordinates including the polar coordinates for 2D systems may be used for 
the velocity space. Afterwards, the continuous Maxwell equilibrium should also be dis- 
cretized. The last step is to adopt an appropriate numerical scheme for the space and time 
discretization. Therefore , we can se e that LB methodology closely resembles the DVM 
problem solving process. iLud (j2000l ) noticed this similarity and stated "the LB equation 
is essentially DVM with finite discrete velocities and fully discretized space and time tied 
to the discrete velocity set". For simulating rarefied gas flows, this similarity is impor- 
tant as we have shown how the LB framework is developed from the Boltzmann-BGK 
equation. 

For both DVM and LB methods, the most critical task is to discretize the velocity 
space. When the Gauss-Hermite quadrature is used in DVM, the discretization of the 
velocity space in these two methods are the same, which may indicate that the LB models 
with sufhciently accurate Gauss-Hermite quadrature can capture the higher-order non- 
equilibrium effects in the rarefied gas flows. This in deed is confirmed by the simulation 
results presented in FiglU which we will discuss in detail in Section 4. 

However, an important advantage of the LB models is the "stream-collision" mecha- 
nism which is mainly inherited from the lattice gas automata. This "stream-collision" 
mechanism makes the LB method easy to understand and simple for computer program- 
ming. Therefore, the "stream-collision" mechanism is an important feature of the LB 
models which distinguishes them from DVM. The coupled time step and physical space 
in the LB models will dramatically reduce the computational cost. In addition, DVM 
relies heavily on mathematical techniques which depend on specific problem, while the 
LB methodology is straightforward and more suitable for developing a generic simulation 
package for engineering design. 

3.3. Lattice Boltzmann equation and linearized BGK equation 

By introducing tl) to denote the unknown perturbed distribution function and assuming 
the fiow is weakly non-equilibrium, / can be approximated by 

/ = + (3.1) 

where 
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which is the global (absolute) equilibrium distribution function. Using the Taylor series 
to expand the local equilibrium distribution function and keeping the terms up to the 
first order, one can obtain the following equation 



l + ^-u+]^{T-l)[e -D) 



(3.3) 



dt ^ ' " Kn 
where we assume the flow is incompressible. Using Eq. p.ip . we can obtain the linearized 
BGK equation: 

^+^-V^ + g-[V^V'-(l + V')l)] = -;^|^- ^■u+'^{T-\){e-D) |. (3.4) 
For lattice Boltzmann models, one can rewrite Eq. (|2.6p as 

/(r, t) « /^(r, I, t) = [1 + ^(r, t)] , (3.5) 

where 

^ 1 

^(r,|,i)-^-a(")(r,i)x("HO- (3-6) 

Substituting Eq. (|3.5p into the Boltzmann-BGK equation and keeping the first- and 
second-order expansions of the equilibrium distribution, one can obtain 



^ + ^ • + g • [V^^ _ (1 + ^) I] = _ I . , (3.7) 



and 

^+|-Vv.+g-[Ve^ -{l + ^)i] = -±- \^^-^.u-\[{^-uf~u^ + {T- l){e - D)] 

(3.8) 

Because w(^) is equal to /o, we can observe the following interesting facts by comparing 
Eqs. (|3.7[ ) and (|3.8p with Eq. (|3.4p . First of all, by keeping the first order Hermite expan- 
sion, the essential LB model equation is the same with that of the isothermal {T — \) 
linearized BGK equation except the body force term. This implies that Lp is indeed equiv- 
alent to -0 though ip is prescribed to include only the finite order terms of the Hermite 
polynomials (cf. Eg. 112. 6\) and Eg. I13.5\) ). Therefore, the LB equation with the first order 
terms should be as good as the linearized BGK equation for isothermal flows. This indi- 
cates that high-order Hermite expansion is not necessary for rarefied gas flows. Secondly, 
with the second order Hermite expansion, there is an extra velocity term ^ [(^ • it)'^ — u^] 
for the LB equation in comparison to the linearized BGK equation. However, for flows 
with low Mach number, this term is a higher-order small quantity which can be ignored. 
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This is the reason why the Hermite expansion order is reported to make neghgible dif- 



ference on the simulation resuhs (see lKim et a/.ll2008f l. In fact, the first order expansion 



is sufficient to obtain the accurate results for isothermal rarefied flows with low speed. 
Furthermore, the LB equation with the second order expansion can in principle describe 
thermal problems since the temperature information is included in Eg. p.Sp . which at 
least has the same capability as the linearized BGK equation, though the BGK kinetic 
model gives wrong Pr number. Thirdly, the treatment of the body force makes the dif- 
ference between the LB model and the linearized BGK equation. It is because that the 
linearized BGK model keeps the full information while the LB model uses the Hermite 
expansion to approximate V^/, i.e. Vj(/3 — (1 + (^) ^. However, for the problem is not far 
from equilibrium state, this difference is not important, which will be confirmed by the 
numerical simulations in Section 4. 

From the above analysis, we can see that the Hermite expansion order does not de- 
termine th e accu racy of LB models for rarefied gas flows as described by Eq.(4.7) in 



Shan et al 



(|2006f ) . The Hermite expansion provides a means to approximate the equilib- 
rium distribution and the body force in the kinetic equation. Therefore, the LB equation, 
similar to the linearized BGK equation, is an approximation of the Boltzmann-BGK 
equation. In contrast to the Grad's moment method, LB models include the information 
of any order moment though it may not be accurate. For instance, with the flrst order 
expansion, the LB model equation is as the same as the isothermal linearized BGK equa- 
tion in the incompressible limit, which will give accurate results for any order velocity 
moment. When the Mach number of flow increases, high-order terms in the Hermite ex- 
pansion become important. Therefore, the order of Hermite expansion is important to 
simulate compressible flows rather than rarefied flows. 

To capture non-equilibrium effects in rarefied flows, the Gauss-Hermite quadrature 
is the key as it determines the discretization accuracy to the model equation. With 
sufRciently high order of the Gauss-Hermite quadrature, LB models can give excellent 
numerical results, e.g. the results presented in Fig.([T|) where 400 discrete velocities are 
used are identical to the DVM solution. Considering the similarity of the LB equation 
and linearized BGK equation, insufficient quadrature order should be responsible for the 
failure on capturing the constitutive relations in the Knudsen layer because the kinetic 
boundary condition have been well accepted in solving the linearized BGK equation. 

In summary, the LB method is essentially a special discrete velocity model, which 
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approximates the Boltzmann-BGK equation with finite discrete velocities and fully dis- 
cretized space and time tied to the discrete velocity set. The capability of LB equation 
is similar to the linearized BGK equation for simulating rarefied gas flows. The Hermite 
expansion order determines the model equation and is important for compressible flows. 
It has no direct effect on the accuracy of capturing high-order non-equilibrium effects. 
Meanwhile, the Gauss-Hermite quadrature as a discretization technique for the velocity 
space directly determines whether the LB models can describe rarefied flows accurately. 



4. Simulations and discussion 

In addition to the above theoretical analysis, we will numerically evaluate the LB 
models. To exclude the boundary condition effect, we choose the standing-shear-wave 
problem as the be nchmark case, which wa s specially designed for assessing the accuracy 



of various models ( Lockerbv fc Reese 



20081 ). It is a shear flow driven by a temporally and 



spatially oscillating body force, which can be written as the following form: 

= Ae">'^ cosOy, (4.1) 

where Fx is the body force in the direction x (which is perpendicular to the y direction) , 
A is the amplitude, and is the wave number and (j) is the frequency. This isothermal 
problem is sufficiently simple because the flow direction is perpendicular to the space 
variation but it is intended to capture the shear-dominated characteristic of microscale 
flows. Furthermore, the distinct advantage is that the boundary is not important here so 
that one can focus on the model itself without the interference from gas molecule/wall 
interactions. With Eq. (|2.2|l . the body force becomes: 

=ie''^*cosy, (4.2) 

where 9 is considered as a measure of the characteristic length, and A = Another 
distinctive advantage for using this benchmark problem is that analytical solutions can 
be obtained for many hydrodynamic models, such as the Navier Stokes equation and the 
regularized 13-moment model (R13). For convenience, the R13 solution is listed as below: 

{288KnH - 5W(l)Kn^ + (520i - 22502i) Kn'^ - 375(t>Kn + 150i) A 

" ~ ^ 288(l3Kn'^ + {510(j)^i - 270i) Kn^ + (7450 - 22503) Kn^ + (37502i - I50i) Kn -f 1500 ' 

(4.3) 
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where u denotes the velocity amplitude. One can refer to I Locker bv fc Reesd (|2008l ) for 
the detail of hydrodynamics models. 

The discrete velocity method of solving the linearized BGK equation has already bee n 
served as a benchmark for the standing shear wave problem bv lLockerbv fc Reesd (|2008[ ). 
where the linearized BGK model Ea. p.4p can be simply written in the scalar form for 
isothermal flows: 



'at 



dy 



9-0 



- (1 + 



1 



{ixUx - 1p) ■ 



(4.4) 



Since the problem is essentially one-dimensional, one can eliminate by multiplying the 
above equation with — = 
resulting equations are: 



above equation with --;=e ^^/^ and -^^^^se respectively. Integrating over ^^x, the 



dGa 

dt 



dGa 



dy 



1 



dGb , ^ dGb p , 1^ 1 / r \ 
-Q^+^y^~FAGa + l) = j^{u.-Gb), 



where the reduced unknown functions Qa and are defined as 



Ga^ 



1 



Gb^ 



1 



/27r 



The macroscopic velocity can be expressed as 

1 



Gbd^y 



(4.5) 
(4.6) 

(4.7) 
(4.8) 

(4.9) 



To solve Eqs. (l4.5p and (|4.6p . the essential task is to choose an appropriate quadrature 
to discretize the velocity space which ranges from —oo to oo. The typical highly accurate 
choice for low speed rarefied gas flows is the Gauss-Hermite quadrature, which is used 
here. Based on the discretization of the phase space, the integration operation over the 
velocity space is converted to sum operation and then a series of equations like Eq. (|2.13p 
are obtained. Naturally, the discretized Maxwell equilibrium distribution can also be 
obtained by directly using its value on the grid of the velocity space. One can then 
use typical numerical methods such as finite difference scheme (e.g., the Lax-Wendroff 
scheme) to solve these equations respectively. 

When the same Gauss-Hermite quadrature with 400 discrete velocities are used in the 
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DVM solution of the linearized BGK equation and our LB model, Fig[T] shows that the 
results for both velocity and shear pressure amplitude are nearly identical for a broad 
range of Knudsen numbers from 0.1 to 1.5. Even with the first order Hermite expansion, 
the LB model can predict shear pressure accurately, which confirms that the Hermite 
expansion order does not directly affect accuracy of the LB models in capturing non- 
equilibrium effects measured by the Knudsen number. 

Although the standard LB model (D2Q9) is not sufficiently accurate in comparison 
with the DVM solution, high-order LB model (D2Q16) with minimal increase of the 
discrete velocity set can produce good results. FiglT] shows that the LB model with 
increasing order of the Gauss-Hermite quadrature can closely approximate the linearized 
BGK equation. Therefore, in comparison with the DVM simulation, LB method can 
provide a practical engineering design simulation tool which can produce reasonably 
accurate results with significantly reduced computational cost. 

As discussed in Section 2.3, at least three factors will influence the problem-solving 
process, i.e. the numerical scheme for solving Ea. (|2.13[) . the order of Hermite expansion 
and Gauss-Hermite quadrature. For the numerical scheme, our second-order scheme is 
essential as discussed in Section 2.2. Regarding the role of Hermite expansion and Gauss- 
Hermite quadrature, we have theoretically proved that the Gauss-Hermite quadrature 
rather than the order of Hermite expansion is key to capturing non-equilibrium effects 
accurately. The numerical simulations have also performed to testify our theoretical anal- 
ysis. 

In Figs [2] and m the simulation results of the three LB models are compared with the 
solutions of directly solving the linearized BGK equation and the Navier Stokes equation. 
The expansion of the equilibrium distribution function and the forcing term is second- 
order for the D2Q9 model, third-order for the D2Q16 and forth-order for D2Q25. The 
results in FigOshow that the prediction for velocity amplitude of the D2Q25 model are 
in excellent agreement with the DVM solution of the linearized BGK equation across a 
broad range of Knudsen number [Kn e [0,1.5]) for the quasi-steady and time- varying 
problems with 6 up to 0.25. Meanwhile, the results of the D2Q9 model deviate from 
the DVM solution of the linearized BGK equation significantly. Surprisingly, the D2Q9 
model does not agree with the results predicted by the Navier Stokes equation. Figl3] 
shows the velocity wave phase lag, which suggests that high-order LB models perform 
better in the transition flow regime. 
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Although FigsO and [3] demonstrate that increasing order of LB model in terms of the 
Hermite expansion and Gauss-Hermite quadrature will lead to more accurate results, we 
still do not know the exact role the orders of the Hermite expansion and the Gauss- 
Hermite quadrature play. Therefore, we single out the effect of the Hermite expansion in 
FigUl where the results of the LB models with the same quadrature but different Hermite 
expansion order are compared. The results clearly show that the Hermite expansion order 
for both the force and the equilibrium distribution function make negligible difference 
to the simulation results. Even the first order expansion is sufficient to obtain accurate 
velocity for the D2Q25 model. The simulation results support our theoretical analysis 
that the Hermite expansion has no direct influence on model accuracy for capturing non- 
equilibrium effects. Specifically, the LB model equation determined by the first order 
Hermite expansion is sufficient for a typical gas flow in micro-devices where the Mach 
number is usually small. In contrast, the Gauss-Hermite quadrature determines the model 
accuracy as the higher-order quadratures give better results. 

Not only the order of quadrature but also the abscissae may influence the model 
accuracy. Therefore, the simulation results of the three LB models with the same order 
quadrature but different abscissae are compared in FigO Although increasing quadrature 
order will lead to improved accuracy, more discrete velocities may not improve the model 
performance if the quadratures are the same order. For example, the quadratures of the 
D2Q16, D2Q17 and D2Q21 models are the same order. Surprisingly, the D2Q16 model 
produces the results better than the other two models with more discrete velocities. 
The reason may be attributed to that the abscissae of the D2Q16 model has better 
symmetry. In addition, all these models are better than the D2Q9 model which has low 
order quadrature. Therefore, appropriate abscissae may improve the model accuracy and 
reduce the computational costs with smaller number of discrete velocities. 



Since iLockerbv fc Reesd (|2008l ) has shown that the R13 equation gives the best perfor- 
mance among the extended hydrodynamic models, we compare the LB models with the 
R13 model here. Fig[6] shows that, in comparison with the data obtained from directly 
solving the linearized BGK equation, the high-order LB models including the D2Q16 
and D2Q25 models can give better results than the R13 equation over a broad range of 
Knudsen numbers. Therefore, the high-order LB models with modest number of discrete 
velocity set, such as the D2Q16 and D2Q25 models, can offer close approximation to the 
linearized BGK equation. Most importantly, these high-order LB models achieve such 



Jianping Meng and Yonghao Zhang 19 

degree of accuracy at a fraction of computational costs associated with directly solving 
the linearized BGK equation. 

5. Concluding remarks 

We have theoretically and numerically analyzed the high-order LB models for rarefied 
gas flows. The lattice Boltzmann equation is shown to be equivalent to the linearized 
BGK equation in the incompressible limit. When the same Gauss-Hermite quadrature 
is used, both LB and DVM simulations produce results in excellent agreement across 
a broad range of the Knudscn numbers. This suggests the importance of the Gauss- 
Hermite quadrature and the great potential of the LB method for modeling rarefied 
gas flows. While the Gauss-Hcrmitc quadrature is of the most importance to capturing 
non-equilibrium effects, the first-order Hermite expansion on the equilibrium distribution 
function is sufficient to obtain the correct moments for isothermal flows e.g. increasing 
the Hermite expansion order further will not improve the model accuracy. For the same 
order Gauss-Hermite quadratures, the chosen abscissae will influence the model accuracy 
and more discrete velocities may not lead to improved model accuracy. 

Overall, we have demonstrated that LB method offers a computationally efhcient ap- 
proach to solve the BGK equation. We can choose a suitable LB model to meet different 
requirement on model accuracy and computational efficiency, which offers an ideal flexi- 
ble engineering design simulation tool to be able to simulate flows in the continuum and 
transition regimes. 
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Figure 1. The results of D2Q400, D2Q16 and D2Q9 models for the quasi-steady standing shear 
wave (a) velocity wave amplitude, (b) shear pressure wave amplitude. The first-order Hermite 
expansion is adopted for the D2Q400 model. Since the Hermite polynomials for the D2Q400 
model give irrational roots, the Lax-Wendroff scheme is used to solve Ea. (|2.13|l here. The results 
show that the LB model with sufficiently large discrete velocity sets can be very accurate. 
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Figure 2. Velocity wave-amplitude as a function of the Knudsen number, where the expansion 
order for the equilibrium distribution function and the force term is N which is 2, 3, 4 for the 
D2Q9, D2Q16 and D2Q25 models respectively, and the order of Gauss-Hermite quadrature is 
2N + 1. 



Jianping Meng and Yonghao Zhang 



25 



<f=0.Z5 



BGK 
D2Q9 
D2Q16 
D2Q25 
Navier Stokes 
R13 




Figure 3. Velocity wave phase lag as a function of the Knudsen number, where the expansion 
order for the equilibrium distribution function and the force term is N which is 2, 3, 4 for the 
D2Q9, D2Q16 and D2Q25 models respectively, and the order of Gauss-Hermite quadrature is 
2JV + 1. 
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Figure 4. Velocity wave-amplitude varying with the Knudsen number. The models are named 
according to the rule D2Qn - Yth where n denotes the number of discrete velocities, Y the 
expansion order for the equilibrium expansion and the force term. 
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Figure 5. Velocity wave-amplitude varying with the Knudsen number, where the three 
models with the same order of quadratures but different abscissae are compared. 
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Figure 6. Velocity wave-amplitude varying with the Knudsen number where the results of LB 
models are compared against the solution of the R13 equation. 



